{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "friendly-bloom",
   "metadata": {},
   "source": [
    "<center><h1>LU分解与列主元的Gauss消去法解线性方程组<h1/><center/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "current-visibility",
   "metadata": {},
   "source": [
    "## 实验内容\n",
    "用LU分解与列主元的Gauss消去法解线性方程组\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\-3.0&2.099999&6.0&2.0\\\\5.0&-1.0&5.0&-1.0\\\\2.0&1.0&0.0&2.0\\end{pmatrix}\n",
    "\\begin{pmatrix}x_1\\\\x_2\\\\x_3\\\\x_4\\end{pmatrix}=\\begin{pmatrix}8.0\\\\5.900001\\\\5.0\\\\1.0\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "(1)用$A=LU$分解，求$L,U$,并求解向量$x$\n",
    "\n",
    "(2)用列主元的Gauss消去法求解向量$x$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "directed-minority",
   "metadata": {},
   "source": [
    "## 实验原理\n",
    "\n",
    "### LU分解\n",
    "\n",
    "直接三角分解的计算公式\n",
    "\n",
    "计算$U$的第$r$行,$L$的第$r$列元素($r=1,2,\\cdots,n$)\n",
    "$$\n",
    "u_{ri}=a_{ri}-\\sum_{k=1}^{r-1}l_{rk}u_{ki}\\qquad i=r,r+1,\\cdots,n;\\\\\n",
    "l_{ir}=(a_{ir}-\\sum_{k=1}^{r-1}l_{ik}u_{kr})/u_{rr}\\qquad i=r,r+1,\\cdots,n,且r\\neq n\n",
    "$$\n",
    "求解$Ly=b,Ux=y$的计算公式：\n",
    "$$\n",
    "\\begin{cases}\n",
    "y_1=b_1,\\\\\n",
    "y_i=b_i-\\sum_{k=1}^{i-1}l_{ik}y_k \\qquad i=2,3,\\cdots,n\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "x_n=y_n/u_{nn},\\\\\n",
    "x_i=(y_i-\\sum_{k=i+1}^n u_{ik}x_k)/u_{ii} \\qquad i=2,3,\\cdots,n\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dense-standing",
   "metadata": {},
   "source": [
    "### 列主元的Gauss消去法\n",
    "\n",
    "1. 对于$i=1,2,\\cdots,n-1$\n",
    "\n",
    "    1. 选取列主元\n",
    "       $$\n",
    "       |a_{i_m,i}|=\\max_{i\\leq k\\leq n}|a_{ki}|\n",
    "       $$\n",
    "\n",
    "    2. 若列主元$a_{i_m,i}=0$，则停止\n",
    "\n",
    "    3. 若$iₘ \\neq i$,则交换行\n",
    "       $$\n",
    "       a_{ij}\\leftrightarrow a_{i_m,j}\\\\\n",
    "       b_i\\leftrightarrow b_{i_m}\n",
    "       $$\n",
    "\n",
    "    4. 消元计算\n",
    "\n",
    "       对于$k=i+1,\\cdots,n$\n",
    "       $$\n",
    "       a_{kj}\\leftarrow a_{kj}-\\frac{a_{ki}}{a_{ii}} a_{ij}\\\\\n",
    "       b_{k}\\leftarrow b_k - \\frac{a_{ki}}{a_{ii}}b_i\n",
    "       $$\n",
    "\n",
    "2. 回代求解\n",
    "    1. $b_n\\leftarrow b_n/a_{nn}$\n",
    "    2. 对于$i=n-1,\\cdots,2,1$\n",
    "        $$\n",
    "        b_i\\leftarrow (b_i-\\sum_{j=i-1}^n a_{ji}b_j)/a_{ii}\n",
    "        $$\n",
    "\n",
    "最终$x\\leftarrow b$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "natural-algeria",
   "metadata": {},
   "source": [
    "## 编程实现、计算实例、数据、结果\n",
    "> 用Julia语言实现"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constant-russell",
   "metadata": {},
   "source": [
    "先导入依赖和输入计算所需的数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "trained-encyclopedia",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra,BenchmarkTools,LambdaFn\n",
    "\n",
    "include(\"../Code/MatrixKit.jl\")\n",
    "include(\"../Code/Swap.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "private-farming",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$A=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\-3.0&2.099999&6.0&2.0\\\\5.0&-1.0&5.0&-1.0\\\\2.0&1.0&0.0&2.0\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$A=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\-3.0&2.099999&6.0&2.0\\\\5.0&-1.0&5.0&-1.0\\\\2.0&1.0&0.0&2.0\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mA=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\-3.0&2.099999&6.0&2.0\\\\5.0&-1.0&5.0&-1.0\\\\2.0&1.0&0.0&2.0\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "architectural-touch",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$b=\\begin{pmatrix}8.0\\\\5.900001\\\\5.0\\\\1.0\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$b=\\begin{pmatrix}8.0\\\\5.900001\\\\5.0\\\\1.0\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mb=\\begin{pmatrix}8.0\\\\5.900001\\\\5.0\\\\1.0\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex b = [8;5.900001;5;1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "temporal-sixth",
   "metadata": {},
   "source": [
    "### LU分解方法计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "second-brunei",
   "metadata": {},
   "outputs": [],
   "source": [
    "function LU(A::Matrix)\n",
    "    n = size(A)[1]\n",
    "    @assert n == size(A)[2] \"The Matrix must be square matrix!\"\n",
    "    L = eye(n)\n",
    "    U = zeros(n,n)\n",
    "    for r = 1:n\n",
    "        U[r:r,r:n] = A[r:r,r:n] - L[r:r,1:r-1]*U[1:r-1,r:n]\n",
    "        L[r+1:n,r:r] = (A[r+1:n,r:r] - L[r+1:n,1:r-1]*U[1:r-1,r:r])/U[r,r]\n",
    "    end\n",
    "    return L,U\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "commercial-expense",
   "metadata": {},
   "outputs": [],
   "source": [
    "function LinearSolveWithLU(L::Matrix,U::Matrix,b::Vector)\n",
    "    n = length(b)\n",
    "    y = zeros(n)\n",
    "    for i = 1:n\n",
    "        y[i:i,1:1] = b[i:i,1:1] - L[i:i,1:i-1]*y[1:i-1,1:1]\n",
    "    end\n",
    "    x = zeros(n)\n",
    "    for i = n:-1:1\n",
    "        x[i:i,1:1] = (y[i:i,1:1]-U[i:i,i+1:n]*x[i+1:n,1])/U[i,i]\n",
    "    end\n",
    "    return x\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "chubby-saint",
   "metadata": {},
   "source": [
    "#### 计算结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "revolutionary-major",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$L=\\begin{pmatrix}1.0&0.0&0.0&0.0\\\\-0.3&1.0&0.0&0.0\\\\0.5&-2.499999999650555e6&1.0&0.0\\\\0.2&-2.3999999996645334e6&0.9599996800001067&1.0\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$L=\\begin{pmatrix}1.0&0.0&0.0&0.0\\\\-0.3&1.0&0.0&0.0\\\\0.5&-2.499999999650555e6&1.0&0.0\\\\0.2&-2.3999999996645334e6&0.9599996800001067&1.0\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mL=\\begin{pmatrix}1.0&0.0&0.0&0.0\\\\-0.3&1.0&0.0&0.0\\\\0.5&-2.499999999650555e6&1.0&0.0\\\\0.2&-2.3999999996645334e6&0.9599996800001067&1.0\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L,U=LU(A);\n",
    "@latex L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "faced-framework",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$U=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\0.0&-1.000000000139778e-6&6.0&2.3\\\\0.0&0.0&1.5000004997903332e7&5.749998499196276e6\\\\0.0&0.0&0.0&5.079998907178727\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$U=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\0.0&-1.000000000139778e-6&6.0&2.3\\\\0.0&0.0&1.5000004997903332e7&5.749998499196276e6\\\\0.0&0.0&0.0&5.079998907178727\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mU=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\0.0&-1.000000000139778e-6&6.0&2.3\\\\0.0&0.0&1.5000004997903332e7&5.749998499196276e6\\\\0.0&0.0&0.0&5.079998907178727\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "figured-basics",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-762.0000900767544"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "detA = (prod∘diag)(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "golden-consolidation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$x=\\begin{pmatrix}-6.09103523174781e-10\\\\-1.0000000008881784\\\\1.0000000000483817\\\\0.9999999998737863\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$x=\\begin{pmatrix}-6.09103523174781e-10\\\\-1.0000000008881784\\\\1.0000000000483817\\\\0.9999999998737863\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mx=\\begin{pmatrix}-6.09103523174781e-10\\\\-1.0000000008881784\\\\1.0000000000483817\\\\0.9999999998737863\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex x=LinearSolveWithLU(L,U,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "maritime-member",
   "metadata": {},
   "source": [
    "#### 性能"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "indie-bacon",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  4.34 KiB\n",
       "  allocs estimate:  46\n",
       "  --------------\n",
       "  minimum time:     2.781 μs (0.00% GC)\n",
       "  median time:      2.908 μs (0.00% GC)\n",
       "  mean time:        3.655 μs (6.24% GC)\n",
       "  maximum time:     317.314 μs (98.53% GC)\n",
       "  --------------\n",
       "  samples:          10000\n",
       "  evals/sample:     9"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark LinearSolveWithLU(L,U,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "valued-finance",
   "metadata": {},
   "source": [
    "#### 误差比较\n",
    "\n",
    "$$\n",
    "\\delta A = LU -A\n",
    "$$\n",
    "\n",
    "其中$L,U$为计算结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "downtown-princess",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$δA=\\begin{pmatrix}0.0&0.0&0.0&0.0\\\\0.0&0.0&0.0&-2.220446049250313e-16\\\\0.0&1.1102230246251565e-16&0.0&0.0\\\\0.0&4.440892098500626e-16&6.4116548942624e-10&0.0\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$δA=\\begin{pmatrix}0.0&0.0&0.0&0.0\\\\0.0&0.0&0.0&-2.220446049250313e-16\\\\0.0&1.1102230246251565e-16&0.0&0.0\\\\0.0&4.440892098500626e-16&6.4116548942624e-10&0.0\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mδA=\\begin{pmatrix}0.0&0.0&0.0&0.0\\\\0.0&0.0&0.0&-2.220446049250313e-16\\\\0.0&1.1102230246251565e-16&0.0&0.0\\\\0.0&4.440892098500626e-16&6.4116548942624e-10&0.0\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex δA = L*U - A"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "color-edinburgh",
   "metadata": {},
   "source": [
    "> 可见误差还是比较小"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "southwest-michigan",
   "metadata": {},
   "source": [
    "$$\\delta b = Ax-b$$\n",
    "\n",
    "其中$x$为计算结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "inclusive-consensus",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$δb=\\begin{pmatrix}0.0\\\\0.0\\\\-1.7892167747390886e-9\\\\-2.3588129227647414e-9\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$δb=\\begin{pmatrix}0.0\\\\0.0\\\\-1.7892167747390886e-9\\\\-2.3588129227647414e-9\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mδb=\\begin{pmatrix}0.0\\\\0.0\\\\-1.7892167747390886e-9\\\\-2.3588129227647414e-9\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex δb = A*x - b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "asian-chemistry",
   "metadata": {},
   "source": [
    "误差$δb$在$10e-9$量级，还是比较小"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stone-arthur",
   "metadata": {},
   "source": [
    "### 列主元的Gauss消去法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fluid-greeting",
   "metadata": {},
   "outputs": [],
   "source": [
    "function GaussEliminate(A::Matrix,b::Vector)\n",
    "    n = length(b)\n",
    "    Ab = [A b]\n",
    "    for i = 1:n-1\n",
    "        iₘ = i-1+argmax(abs.(Ab[i:n,i]))\n",
    "        iₘ == i || @swap Ab[i,i:end],Ab[iₘ,i:end]\n",
    "        Ab[i+1:n,i+1:end] -= Ab[i+1:n,i:i]*Ab[i:i,i+1:end]/Ab[i,i]\n",
    "    end\n",
    "    x = Ab[:,n+1]\n",
    "    for i = n:-1:1\n",
    "        x[i] /= Ab[i,i]\n",
    "        x[1:i-1] -= Ab[1:i-1,i]*x[i]\n",
    "    end\n",
    "    return x\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "disturbed-richmond",
   "metadata": {},
   "source": [
    "#### 计算结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "original-match",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$x=\\begin{pmatrix}2.6645352591003756e-16\\\\-0.9999999999999997\\\\0.9999999999999999\\\\1.0\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$x=\\begin{pmatrix}2.6645352591003756e-16\\\\-0.9999999999999997\\\\0.9999999999999999\\\\1.0\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mx=\\begin{pmatrix}2.6645352591003756e-16\\\\-0.9999999999999997\\\\0.9999999999999999\\\\1.0\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex x = GaussEliminate(A,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "extended-protest",
   "metadata": {},
   "source": [
    "#### 性能"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "hispanic-evans",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  4.86 KiB\n",
       "  allocs estimate:  44\n",
       "  --------------\n",
       "  minimum time:     2.718 μs (0.00% GC)\n",
       "  median time:      2.844 μs (0.00% GC)\n",
       "  mean time:        3.515 μs (8.76% GC)\n",
       "  maximum time:     466.047 μs (98.40% GC)\n",
       "  --------------\n",
       "  samples:          10000\n",
       "  evals/sample:     9"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark GaussEliminate(A,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "elect-warrant",
   "metadata": {},
   "source": [
    "#### 误差比较\n",
    "\n",
    "$$\\delta b = Ax-b$$\n",
    "\n",
    "其中$x$为计算结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "failing-seller",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$δb=\\begin{pmatrix}0.0\\\\0.0\\\\8.881784197001252e-16\\\\8.881784197001252e-16\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$δb=\\begin{pmatrix}0.0\\\\0.0\\\\8.881784197001252e-16\\\\8.881784197001252e-16\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mδb=\\begin{pmatrix}0.0\\\\0.0\\\\8.881784197001252e-16\\\\8.881784197001252e-16\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex δb = A*x - b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faced-johnston",
   "metadata": {},
   "source": [
    "误差$δb$在$10e-16$量级，还是比较小"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "compliant-prison",
   "metadata": {},
   "source": [
    "### *列主元的LU分解*\n",
    "> 补充内容"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "vital-stage",
   "metadata": {},
   "outputs": [],
   "source": [
    "function P⁻¹LU(A::Matrix)\n",
    "    # deepcopy otherwise you will change the origin data\n",
    "    A = A |> deepcopy\n",
    "    n = size(A)[1]\n",
    "    @assert n == size(A)[2] \"The Matrix must be square matrix!\"\n",
    "    Iₚ = zeros(Int8,n)\n",
    "    for r = 1:n\n",
    "        #1 calculate s\n",
    "        A[r:n,r:r] -= A[r:n,1:r-1]*A[1:r-1,r:r]\n",
    "        \n",
    "        #2 choose main element\n",
    "        iᵣ = r-1+argmax(abs.(A[r:n,r]))\n",
    "        Iₚ[r] = iᵣ\n",
    "        \n",
    "        #3 swap row of A\n",
    "        r == iᵣ || @swap A[r,:],A[iᵣ,:]\n",
    "        \n",
    "        #4 calculate L,U same as function LU\n",
    "        A[r+1:n,r:r] /= A[r,r]\n",
    "        A[r:r,r+1:n] -= A[r:r,1:r-1]*A[1:r-1,r+1:n]\n",
    "    end\n",
    "    return A,Iₚ\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "intended-partition",
   "metadata": {},
   "outputs": [],
   "source": [
    "function LinearSolveWithP⁻¹LU(LUᵤₙᵢₒₙ ::Matrix,Iₚ::Vector,b::Vector)\n",
    "    b = b |> copy\n",
    "    n = length(Iₚ)\n",
    "    for i = 1:n-1\n",
    "        i == Iₚ[i] || @swap b[i],b[Iₚ[i]]\n",
    "    end\n",
    "    for i = 2:n\n",
    "        b[i] -= LUᵤₙᵢₒₙ[i,1:i-1]'*b[1:i-1]\n",
    "    end\n",
    "    b[n] /= LUᵤₙᵢₒₙ[n,n]\n",
    "    for i = n-1:-1:1\n",
    "        b[i] -= LUᵤₙᵢₒₙ[i,i+1:n]'*b[i+1:n]\n",
    "        b[i] /= LUᵤₙᵢₒₙ[i,i]\n",
    "    end\n",
    "    return b\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "novel-fantasy",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$LUᵤₙᵢₒₙ=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\0.5&2.5&5.0&-1.5\\\\-0.3&-4.000000000559112e-7&6.000002&2.2999994\\\\0.2&0.9600000000000002&-0.7999997333334223&5.079998906667031\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$LUᵤₙᵢₒₙ=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\0.5&2.5&5.0&-1.5\\\\-0.3&-4.000000000559112e-7&6.000002&2.2999994\\\\0.2&0.9600000000000002&-0.7999997333334223&5.079998906667031\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mLUᵤₙᵢₒₙ=\\begin{pmatrix}10.0&-7.0&0.0&1.0\\\\0.5&2.5&5.0&-1.5\\\\-0.3&-4.000000000559112e-7&6.000002&2.2999994\\\\0.2&0.9600000000000002&-0.7999997333334223&5.079998906667031\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LUᵤₙᵢₒₙ,Iₚ=P⁻¹LU(A);\n",
    "@latex LUᵤₙᵢₒₙ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "circular-picking",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$Iₚ=\\begin{pmatrix}1&3&3&4\\end{pmatrix}$$^T$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$Iₚ=\\begin{pmatrix}1&3&3&4\\end{pmatrix}$$^T$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mIₚ=\\begin{pmatrix}1&3&3&4\\end{pmatrix}\u001b[39m\u001b[35m^T\u001b[39m</center>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex Iₚ T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "revolutionary-alberta",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "<center>$x=\\begin{pmatrix}2.6645352591003756e-16\\\\-0.9999999999999997\\\\0.9999999999999999\\\\1.0\\end{pmatrix}$</center>\n",
       "\n"
      ],
      "text/markdown": [
       "<center>$x=\\begin{pmatrix}2.6645352591003756e-16\\\\-0.9999999999999997\\\\0.9999999999999999\\\\1.0\\end{pmatrix}$</center>\n"
      ],
      "text/plain": [
       "  <center>\u001b[35mx=\\begin{pmatrix}2.6645352591003756e-16\\\\-0.9999999999999997\\\\0.9999999999999999\\\\1.0\\end{pmatrix}\u001b[39m</center>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@latex x = LinearSolveWithP⁻¹LU(LUᵤₙᵢₒₙ,Iₚ,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "proof-landscape",
   "metadata": {},
   "source": [
    "## 结果分析与比较\n",
    "根据上面的计算结果\n",
    "\n",
    "|          | 直接$LU$分解法 | 列主元的$Gauss$消去法 |\n",
    "| :--------: | :--------------: | :---------------------: |\n",
    "| 误差$δb$ | $10e{-9}$量级  | $10e{-16}$量级        |\n",
    "| 性能     | 2.908 μs       | 2.655 μs              |\n",
    "\n",
    "> 这里性能采用运行`中位数时间(median time)`度量\n",
    "\n",
    "### 分析结论:\n",
    "\n",
    "计算精确度：列主元的$Gauss$消去法比直接$LU$分解更高\n",
    "\n",
    "性能：列主元的$Gauss$消去法比直接$LU$分解更好\n",
    "\n",
    "> 这与理论分析的结论一致"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exciting-permission",
   "metadata": {},
   "source": [
    "## 参考文献\n",
    "- [1] 李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社,2018.12\n",
    "- [2] Julia.Docs:Linear Algebra[EB/OL].https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "reasonable-harris",
   "metadata": {},
   "source": [
    "## Copyright\n",
    "\n",
    "> Copyright 2021 by Algebra-FUN(樊一飞)\n",
    ">\n",
    "> ALL RIGHTS RESERVED."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,jl"
  },
  "kernelspec": {
   "display_name": "Julia",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.0"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
